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Abstract. The variational free-Lagrange (VFL) method for shallow water 
is a free-Lagrange method with the additional property that it preserves the 
variational str ucture of sha l low w ater. The VFL method was first derived in 
this context by[Augcnbaum ( 1984|) who discretized Hamilton's action principle 
with a free-Lagrange data structure. The purpose of this article is to assess the 
long-time conservation properties of the VFL method for regularized shallow 
water which are useful for climate simulation. Long-time regularized shallow 
water simulations show that the VFL method exhibits no secular drift in the 
(i) energy error through the application of symplectic integrators; and (ii) the 
potential vorticity error through the construction of discrete curl, divergence 
and gradient operators which satisfy semi-discrete divergence and potential 
vorticity conservation laws. These diagnostic semi-discrete equations augment 
the description of the VFL method by characterizing the evolution of its re- 
spective irrotational and solenoidal components in the Lagrangian frame. Like 
the continuum equations, the former exhibits a div 2 U term which indicates 
that the flow has a very strong tendency towards a purely rotational state. 

Numerical results show (i) the preservation of shape and strength of an 
initially radially symmetric vortex pair in purely rotational regularized shallow 
water and (ii) how the Voronoi diagram retains the history of the flow field 
and (iii) that energy is conserved to C(A 2 ) and potential vorticity error to 
within 5% with no secular growth over a 50 year period. 



1. Introduction 

Motivation. The derivation and analysis of conservative numerical methods for 
global climate modeling is a challenging and active area of research. The traditional 
approach to deriving such methods requires the numerical analyst or engineer to ex- 
ercise some level of their own bias in choosing which conservation laws to enforce in 
the discrete model. To remedy this and systema tically derive th e numerical method 
from a conferred mathematical understanding, ISalmonl (Il983l) pursues the varia- 



tional approach for geophysical modeling. The power of the variational approach 
rests upon Noether's theorem which states that each continuous symmetry of Hamil- 
ton's action principle exhibits a corresponding conservation law. Consequently, the 
numerical analyst need only choose a suitable discretization of Hamilton's action 
principle, with the confidence that not only will the resulting semi-discrete equa- 
tions of motion preserve geometric structure, manifesting in energy conservation, 
but that these equations will conserve additional quantities corresponding to the 
symmetries of the discrete action principle. 
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The free-Lagrange method. lAugenbauml <| 1984T ) presents a variational extension of 
the free-Lagrange method for rotating shallow water- a method which by itself ex- 
hibits a number of favorable properties. In contrast with the majority of methods 
for shallow water, which use a fixed grid for transport computations, the free- 
Lagrange method builds the transport into the evolution of particle positions. This 
step avoids introducing diffusion into the transport- an artefact of mixing the trans- 
port computed from fluxes between cells with the remaining fluid. Consequently, 
the free-Lagrange method exhibits local transport conservation laws such as fluid 
mass and tracer concentration. 

Secondly, the free-Lagrange method represents the layer thickness over a Voronoi 
diagr am. This appro a ch is attractive from both a theoretical and practical perspec- 
tive. ISerrano et all (|2005t ) shows that use of a Voronoi diagram in Hamilton's 
discrete action principle permits translational and rotational symmetries and hence 
respective linear and angular momentum conservation laws for the discrete Eu- 
ler equations. The Voronoi diagram avoids mesh-tangling problems which seriously 
compromise the performance of fixed connectivity mesh based Lagrangian methods. 
This does not mean that use of a Voronoi diagram is not without its own compu- 
tational challenges. The generation of an entirely new Voronoi diagram at each 
time step is computationally prohibitive in three dimensions. Also, computational 
complexity manifests from the absence of constraints on the mesh connectivity. 

Some of these computational efficiency shortcomings hav e been partia l ly ad - 



dressed by alternative approaches cited in the literature. lHarlen et al.l (|1995l ) 



present a split Eulerian-Lagrangian scheme for viscoelastic fluids that retains the 
nodes as material points and reconnects them to produce the optimal Delaunay 
triangulation. This offers an advantage over the free-Lagrange method in that the 
discrete fl uid equations can then be solved using Galerkin finite element approx- 



imations. IPetera and Nassehil (|1996l ) presents a Galerkin-Lagrange finite element 
approximation of a shallow water model for tidal flow which attempt to resolve 
fluctuations between dry and wet regions. Their basis for node adjustment is a 
physical one. The nodes of the mesh are only moved if the neighbourhood of the 
node changes state between wetted and dried. This approach avoids excessive mesh 
distortions but requires a smoothing operator to eliminate numerical oscillations in- 
troduced by succesive modifications of the mesh in very shallow regions. 
Variational methods. Despite their computational efficiency, neither of these two 
previously described approaches preserve variational structure- a distinguish i ng fea - 
ture of the variational free-Lagrange (VFL) method derived by I Augenbaum (1984). 
We trac e the devel o pmen t of variational numerical methods in geophysical modeling 
back to iBunemanl (|1982j ). who describes a two-dimensional Eulerian model based 
upon the Clebsch representation of Hamilton's principle. He advocates a method 
whose stability appears to be contingent on the invariance of the discrete Hamilton's 
principle under particle relabelling. It is now well known, however, that the discrete 
Lagrangian is not invariant under conti nuous transformations of t he labels - there 
is no particle relabelling symmetry (see iRipal . 1 1 9 8 it ISalmonl . Il982t ) . The symmetry 
gives Ertel's theorem of conservation of potential vorticity which in turn recovers 
the known connection with Kelvin's circulation theorem along surfaces of constant 
entropy. In the pres e ntatio n of a numerical analogue of a variational shallow water 
blob model, Salmon ( 19831) conjectures that potential vorticity conservation might 
still be conserved through modifying the discrete variational principle. Specifically 
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he suggests constraining the angular component of the particle velocities. This ap- 
proach is difficult to efficiently implement, however, and remains unsubstaniated. 

More computationally feasible variational approaches have emerged from clas- 
sical Lagrangian numerical methods. Most notably, variational form ulations of 



smoo th particle hydrodynamics (SPH) methods have been developed (see lBonet and Rodriguez-Paz , 
2005, and the references by the same first author therein). In this work, the authors 
actually make use of the variational principle to extend the classical SPH method 
by deriving spatially dependent smoothing lengths (a limiting feature in conven- 
tional SPH methods is the use of a fixed smoothing length). There appears to be 
a notable advantage in using a variational SPH for highly compressible fluids. 
Symplectic integrators. The literature that we have surveyed on the VFL method 
does not consider preservation of the symplectic structure under time discretiza- 
tion of the particle Euler-Lagrange equations. Our primary contribution is to ex- 
ploit the symplectic structure preserved by the dispersively regularized shallow 
water equations. Disp ersive regularizat i on avo ids loss of symplecticity of symplec- 
tic time-steppers (see iDixon and Reichl [2004, for a review of symplectic methods 
for particle methods) by slowing high frequency error modes which would other- 
wise cause numerical instabilities. We demonstrate by numerical experiment that 
the combination of the regularized VFL method with an explicit symplectic time- 
stepper exhibits no secular drift in the energy error for moderate time steps and 
coarse spatial resolutions. This property is desirable for long-time geophysical sim- 
ulations which rely on moderate time step sizes and coarse spatial resolutions for 
computational tractability. 

Overview. We shall begin by revisiting the free-Lagrange method and subsequently 
deriving the semi-discrete VFL equations for regularized shallow water in Section 
[3] In Section [4] we present the main feature of the paper, an explicit symplectic in- 
tegrator for preserving the symplectic structure of the semi-discrete VFL equations 
for regularized shallow water. The form of the regularization operator and its im- 
plementation is described in Section [5j Using the semi-discrete VFL equations, we 
derive the semi-discrete vorticity and divergence equations in Sections 16.11 and 16.21 
and show that the potential vorticity is only conserved if the discrete curl operator 
is judiciously chosen so that its operation on the layer thickness gradient vanishes. 
In Section [3 we present numerical results which show the scaling laws of potential 
vorticity and energy. The VFL method for ID rotating shallow water is shown to 
conserve energy to the order of the integrator and exhibit the geostrophic adjust- 
ment mechanism of rotating shallow water. Numerical simulations illustrate the 
motion of a vortex pair in 2D rotating shallow water and the configuration of the 
Voronoi diagram. We also show the corresponding energy and potential vorticity 
conservation properties over a 50 year period. Section [8] concludes. 



2. The Free-Lagrange Method 

Consider a particle representation of a fluid at time t in which N 2 particles or 
sites X(t) = {Xi(t), . . . ,Xjv2(t)} are individually labelled by a. Each particle is 
inside a polygon of area A a (t) containing the set of distinct points X^(t) closer 
to X Q (i) than to any other site. The polygon is referred to as a Voronoi cell and 
the set of all these closed cells is the Voronoi diagram. A hexagonal Voronoi cell 
is shown in Figure 1 together with the specification of a local index for referring 
to neighboring particles. Note that the assumption that the sites are distinct does 
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not prohibit the formation of shocks - charac teristics may stil l collide as they are 
not the particle trajectories (see, for example, Whithaml . 1974L for an explanation 
of shock formation). 




FIGURE 1 . A hexagonal (n e — 6) Voronoi cell (outlined with the solid line) 
containing the particle with label a. Each cell edge is of length Az£* and indexed 
by The surrounding Dclaunay triangulation, dual to the Voronoi diagram, is 
shown with the dotted line. 



The free-Lagrange method simply uses the particles to represent the material 
velocity field U Q and the Voronoi diagram to construct the piecewise constant 
cellwise layer thickness ft.™ := /i(X a ,t") which, in discrete time, is computed from 
the cell mass conservation law 



(1) 
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3. A Variational Free-Lagrange method for 2D Shallow Water 

Discretization of Hamilton's action principle for rotating shallow water with this 
free-Lagrange fluid representation is the critical step in deriving a VFL method. 
The semi-discrete material description of Hamilton's action principle for 2D rotating 
shallow water in a f-plane is expressed in terms of the particle material velocities 
and piecewise constant cellwise layer thickness. 

Definition 3.0.1 (The Semi-Discrete Hamilton's Action for Regularized Shallow 
Water). The semi-discrete action principle for rotating shallow water over bottom 
topography in a f-plane is 

(2) 

Sd = \ f dtY,™ a (\U a (t)\ 2 + 2R a -U a ) -V, V = ^^2m a (h(^ a ,t) + 2b, 
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with fixed end points SX a (t a ) = SX a (tb) — 0, Va. h is a dispersive regularization 
of the layer thicknes^. 

Under the f-plane approximation, V x R = /ok, k is the unit vector in the 
opposite direction of gravity and /o is the Coriolis parameter which is given by 
/o = 2|cj|, where u) is the angular velocity of rotation relative to an inertial frame. 
g is the gravitation constant, b a is the bottom topography, assumed to be piecewise 
constant over cell a and V is the potential energy. 

Stationarity of the discrete action principle is equivalent to the semi-discrete 
Euler-Lagrange particle equations 

(3) tl a = g [h] * dni* - / k x U Q , ± a = U a , 

i 

where dn^ :— Al^ is the outward normal vector at edge Pi of cell a, Al^ is the 
edge length and the averaging operator = j^- (m a (- a ) + ti^-^J) evaluates 
a cellwise scalar quantity over this edge as the mass weighted mean of that quantity 
over cell a and cell Pi. These weights are constants since, by construction, the mass 
is conserved in each cell. The right hand side of (j3|) gives the form of the cellwise 
discrete gradient operator 

(4) grad(MX Q ,i)) := - £ dn f% 

i 

which is a discrete analogue of the weak form of the continuum gradient operator 
uniquely defined over the interior of a closed cell whose area tends to zero. This dis- 
crete gradient operator is given in weak form, (|3]) applies pointwise to the material 
velocity vector field which is constructed at particle positions only. 

Remark 3.0.2. We note that the discrete gradient operator ([4]) does not define a 
gradient along the cell edges. The prognostic discrete divergence and potential vor- 
ticity conservation laws corresponding to the weak form of continuum conservation 
laws, given in Section [5[ require a definition of the gradient of the regularized layer 
depth along cell edges. We look instead to the linear interpolant of the layer depth 
constructed over the Delaunay triangulation, dual to the Voronoi diagram, to evalu- 
ate a gradient here for prognostic purposes only. We emphasize that the evaluation 
of the gradient along cell edges just enables a more complete characterization of the 
conservative properties of the VFL method and is not intrinsic to the VFL method. 



4. A Symplectic Time Stepping Scheme 

The semi-discrete Euler-Lagrange particle equations ((3]) are symplectic and con- 
serve energy. A symplectic time-stepper, such as the explicit second order Stormer- 
Verlet scheme, in particle position and canonical momentum co-ordinates (X™ , = 



The details of the dispersive regularization operator are given in Section [5] 
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to q U") at time i" 



P" + *=P£-^x(XS), 



pTl-1 

(5) X" +1 = X" + At— 



p^ = pS +i -^(xs +1 ), 

preserves the symplectic two-form uj n = Yl a =i C ^-S ^ rfP„, where Vx denotes the 
gradient of the potential energy. 

Discussion. Integrators which conserve the Hamiltonian are often regarded as a 



strong manifestation of unconditional numerical stability ([Lewis and Simol . I1994T ) 
and therefore permit the use of large time steps. Symplectic integrators, on the 
other hand, conserve a modified Hamiltonian which is derived by backward error 
analysis. When used as time steppers for variational free-Lagrange methods, both 
the particle resolution and the time step size must be chosen to ensure proximity 
of the modified Hamiltonian to the exact Hamiltonian. Otherwise, particles collide 
resulting in singularities in the Jacobian of the surrounding Voronoi cells. In such 
cases, the particles must be remapped and the variational structure is lost. 

Dispersive regularization of the layer depth prevents particles colliding with a 
time step size generally no smaller than required by semi-implicit semi-Lagrangian 
methods. These methods, despite being dissipative, are preferred by practitioners 
as they permit larger time steps than Eulerian methods. The reader is referred 



y p ermit 

to iReic bl (|2006l) for details of dispersive regularization methods for particle mesh 



approximations of shallow water and their comparison to the dispersive properties 
of semi-implicit time stepping. 

The principal advantage of using the variational free-Lagrange method combined 
with symplectic integrators over non-symplectic time steppers for semi-Lagrangian 
methods is the general absence of secular drift in the energy error. The ever diver- 
sifying set of computational models which crucially rely on the ability to tractably 
simulate long-time dynamics without secular numerical dissipation motivates the 
development of new computationally competitive geometric numerical methods 
such as the VFL method. Global climate models are arguably one of the strongest 
candidates in this regard- implementers are challenged by the need to tractably 
resolve complex geophysical dynamics over long periods. 



5. Dispersive Regularization 

The modification of the semi-discrete action principle, to include a dispersively 
regularized potential, ensures that the corresponding time stepping scheme {H} re- 
mains stable for moderately large time steps at computationally tractable resolu- 
tions for long-time global climate simulations. This is frequently accomplished by 
approximating the inverse of a Hclmholtz elliptic operator 



(6) 



h=(l-a 2 V)- 1 h, 
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which includes a characteristic dispersive regularization length scale a. We con- 
struct a conforming finite element approximation of this operator on the Delaunay 
triangulation Td which is dual to the Voronoi diagram over the domain Q. 

More precisely, find hd € V(hj) with a(hd, <j>) = b(hd, </>) for all 4> in the space 
of all linear functions. a(, ) and &(, ) are bilinear forms respectively defined on the 
finite dimensional vector spaces V(hd) and V(hd) where 

(7) V(h d ) := {h d e C(ft) | h\ K g VxQC) V/C e T d c 

Let p^, px. G Pi be the polynomial interpolating functions on triangle K, where 

(8) pic(x,y) =^2h(x,y)Ni(x,y) and p K (x, y) = ^ M x > v) N i( x , v), 

i i 

and the locally indexed basis functions over triangle K, take the form 
Ni{x,y) = [(x 2 y 3 - x 3 y 2 ) + (2/2 - 1fy)x + (x 3 - x 2 )y)]/2A, 

(9) N 2 {x, y) = [(x 3 yi - ssijfe) + (j/ 3 ~ l/i)a; + (a?i - x 3 )y)]/2yl, 
^3(1, y) = [(xm - x 2 yi) + (yi - y 2 )x + (x 2 - x%)y)]/2A, 

where A denotes the area of triangle /C. Defining hd(x,y) = pic{x,y), hd{x,y) — 
pic(x,y) V.x G JC G Td and constructing the sparse stiffness matrix A and mass 
matrix M representing a(hd, Ni) — b(hd, Ni), Vi to give the sparse linear system 

(10) Ah = Mh, 

where the matrices are defined from the weak form of © 

A v = a(N h Nj) - [ d2xN i N o + ® 2 ( N i),k( N i),k, 

(ID 

M l3 = biN^Nj) - / ^xN.Nj. 
K.er d Jk 

Computational procedure. At each time step, the linear system ([T0|) is solved by 
inverting the sparse stiffness matrix to give the regularized layer thickness at each 
Voronoi site. This layer thickness is a prognostic variable used to evaluate the 
discrete cellwise gradient operator (j4]) in the right hand side of (j3|). We outline the 
basic mechanism for applying the regularization to the VFL method: 



(1) The unregularized layer thickness is updated at each time step according 
to the mass conservation law ([1]); 

(2) The linear interpolant of the updated unregularized layer thickness h d is 
constructed over the Delaunay triangulation; 

(3) The Helmholtz linear system (fTU)) is solved for the linear interpolant of the 
regularized layer thickness; and 

(4) The discrete gradient Q is evaluated from the nodal values of the regular- 
ized layer thickness. 



s 
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The linear interpolant of the layer thickness over the triangulation uniquely de- 
fines the gradient of the regularized layer thickness at the edge between neighboring 
sites with global indices a and /3j 



(12) (grad (ho))^ ■= ^ (hp, - h a ) , 

where is the Euclidean distance between sites a and j3{. This form of the 
gradient operator will now be used to derive the diagnostic semi-discrete potential 
vorticity and divergence equations. 



6. The Shallow Water Potential Vorticity and Divergence 

Equations 

Geophysical flows are generally dominated by potential vorticity dynamics, a 
property which together with the divergence underlies the prognostic momentum 
equation. We seek conditions under which the discrete momentum equations have 
underlying potential vorticity and divergence dynamics which are consistent with 
Ertel's theory of potential vorticity conservation and the continuum divergence 
equation respectively. To this end we derive the semi-discrete potential vorticity and 
divergence equations from the discrete momentum equations and identify spurious 
terms arising from the semi-discretization. 

The curl and divergence of velocity respectively kinematically equate to the 
solenoidal and irrotational components of the Hclmholtz decomposed velocity field. 
For this reason, the curl and divergence evolution equations are typically derived to- 
gether and describe how each corresponding component of the velocity field evolves. 
Moreover, the decomposition has the added benefit of largely separating the tem- 
porally fast modes characteristic of the irrotational flow from the temporally slow 
modes characteristic of the solenoidal flow. 

The formulation of the semi-discrete potential vorticity and divergence equa- 
tions depends on the definition of discrete curl, div and grad operators. We have 
already shown that semi-discretization of Hamilton's action principle with the free- 
Lagrange data structure defines the discrete grad operator over the interior of the 
Voronoi cells. The remaining discrete operators, including the definition of the layer 
thickness gr adient along cell edges is not defined by the VFL method. 

Following Ringl er and Randall (|2002h . we define the normal component of the 



discrete curl operator by analogy with the derivation of the discrete divergence 
operator given by e (fT4|) . The Voronoi diagram is assumed to rest in the X-Y plane 
so that the vertical unit vector k is normal to the plane. From Stoke's theorem we 
arrive at 



(13) J- / k.VxUd^ll f-lM = -J-f%ft- / IM, 

A a J A a A a JdAc A a T~ J J 8A% 4 

which approximates to 

( 14 ) T £ ^ • / * A~ E dr « ' : = fe ' curl (u a ), 
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where dr^ ■— AZ^\ The discrete curl operator applies to a cell edgewise vector 
field and not to the vector field constructed at particle positions. As previously 
mentioned, because the discrete gradient operator Q is not uniquely defined over 
cell edges, we instead evaluate the gradient of the layer thickness at the cell edges 
using the interpolant constructed over the dual mesh. 

The additional discrete operators needed to establish discrete divergence and 
potential vorticity laws are the discrete div and curl operators which respectively 
take the form 



(15) 



k.curl(. a ) :=— 5>r*.(.£ 



2 — 1 



and are defined in cell a by evaluating the argument along the cell edges. 

6.1. The semi-discrete potential vorticity equation. The vertical component 
of the discrete curl of the semi-discrete Euler-Lagrange equation over cell a takes 
the form 

(16) k • curl(U Q ) = -gk ■ curl (gr&d(h(X a , t))) - / k • curl(k x U Q ). 

Using the definition of curl, after several calculations provided in Section IA.11 the 
left hand side of (TTBl becomes 



(17) k • curl(U Q ) = ^ (k • dirl(U a )) + div(U«)k • curl(U Q ), 
where 

(18) lTt dT « = Wt ^ Al ^ = A lTt ifl) * = Au * := u * +1 

Substituting the definition of the relative vorticity £ Q = k • curl(U Q ) into the last 
line of (|17l) . we arrive at 

(19) k ■ curl^ = -^( a + Cadiv(U Q ). 

Taking this expression and substituting it into (|16p gives the semi-discrete vorticity 
equation 

(20) ^- f ( a - C Q div(U Q ) + / k • curl (U^) = 0, 

where the right-hand side vanishes because the discrete curl of the edgewise gradient 
(IT2l) evaluates to 



k • curl(grad(/i Q ,)) = ^ (grad(h a )Y ■ dr^ 

a i 

( 21 ) 1 fr 7 \ . ft 



= 0. 
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Using the following properties of the discrete curl and div operators 

, n c , Tie 

k.curl(Ui) - ^-E d ^ ' ( U «) ■= ^E d ^ ' ( U « x £ ) 

( 22 ) Q ^ " 

= i-E dn « • (U Q ) = div(U Q ), 

a i=l 

the reader can easily confirm that (I2P1) simplifies to 

(23) ^ (C„ + /o) + (Co + /o) div(U Q ) = 0. 

Potential vorticity. The semi-discrete potential vorticity equation follows immedi- 
ately from the semi-discrete vorticity equation and the semi-discrete mass continuity 
equation 

(24) ^-h a = -h a diy(V a ). 

Substituting this expression for the discrete divergence of the particle velocity U Q 
into (|23| gives the semi-discrete potential vorticity conservation law 

(25) R (C + f ) + = ± (^) = 0. 
V ' Dt J "' \ h a J Dt Dt\ h a ) 

We observe that this semi-discrete potential vorticity conservation law is condi- 
tional upon the discrete curl and gradient operators satisfying (|2~T|) . 



6.2. The semi-discrete divergence equations. We now repeat the previous 
steps for deriving the semi-discrete potential vorticity equation, only this time, 
taking the discrete divergence of the semi-discrete Euler-Lagrange equations. The 
discrete divergence of the semi-discrete Euler-Lagrange equation over cell a takes 
the form 

(26) div(U Q ) = - 5 div (grad(h(X Q , t))) - / div(k x U Q ). 

Using the definition of div, after several calculations provided in IA.21 the left hand 
side of (l26l) becomes 



(27) 



div(U Q ) = J- J2 ■ dn^ = ^div(U Q ) + div(U Q ) 2 - r a . 



Denoting the divergence of U a as S a :— div(U Q ), (|27p becomes 

(28) div(-^U a ) = ^ a + ^-r a . 
So the semi-discrete divergence equation over cell a is 

(29) ^S a +S 2 a -T a + 5 div(grad(/i Q )) - /odiv(U^) = 0. 

Remark 6.2.1. We compare the semi-discrete shallow water divergence term in 
(f2"5| with its continuous form given by 
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(30) 

Dt 



div ( -^-U ] = d t 8 + div{\J ■ VU) 



= d t 6 + d x (Ud x U + Vd Y U) + d Y {Ud x V + Vd Y V) 

= d t S + {d x Uf + (d Y V) 2 + 2d Y Ud x V + Ud x U + Vd Y V + Vd XY U + Ud XY V 
= d t S + (d x U) 2 + {d Y Vf + 2d Y Ud x V + U ■ WS 

= —S + 5 2 + 2{d x Vd Y U - dxUdyV) 

= ^ + 5 2 - 2de£(VU). 

Intruiguingly, because the Lagrangian form of the continuum shallow water di- 
vergence equation is less visited by the literature, we find the need to emphasize 
the physical significance of the div 2 JJ term, which has a very strong tendency to 
force the fluid back to a purely rotational state. The discrete variant of this term 
also appears in the semi- discrete divergence equations. Through appealing to fur- 
ther analysis of this term 's effect on continuum vortex dynamics can we formulate 
an analogous statement in support of correctly resolved vortex dynamics in simula- 
tion. We also note that T a only evaluates to -j- §„ A U-d(U _L ) = 2c?et(VU) in the 
continuum limit. So like the semi- discrete potential vorticity, we are able to iden- 
tify analogous physically significant terms, but the extent to which these quantities 
are conserved in simulation is governed by the consistency of the choosen discrete 
operators with the approximation solution fields. 



7. Numerical Experiments 

This Section describes three numerical experiments using the VFL to simulate 
the rotating regularized shallow water equations on a periodic domain. The pur- 
pose of the first experiment is to investigate the conservative properties of the 
VFL method for rotating ID shallow water. In the second experiment, we con- 
sider the problem of geostrophic adjustment of rotating ID shallow water and 
show that the VFL method produces result s which are c onsistent with the the- 
ory of geostrophic adjustment established bv lRossbv (1938). This theory describes 



the physical mechanism by which perturbed rotating shallow water recovers to a 
geostrophically balanced state. Layer depth perturbations h' of the scale L' << Ld, 
where Ld = \/gHo/ fo is the Rossby deformation radius for shallow water in a f- 
plane, result in gravity waves which propagate energy and momentum away from 
the source leaving behind a geostrophically balanced flow. It is precisely because 
this mechanism is consistent with the advection law of potential vorticity that geo- 
physical modelers use this experiment to assess the integrity of the discrete shallow 
water velocity and layer depth equations. These two experiments are initialized by 
perturbing the Voronoi cell masses with a Gaussian perturbation and the velocity 
field is initially zero. 

In the third experiment, we simulate the motion of a vortex pair in purely ro- 
tational 2D shallow water flow and monitor the energy and potential vorticity 
conservative properties over a much longer period of 50 years. We verify that the 
geostrophic balances state is preserved - the vortex pair undergoes pure rotation 
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without shape deformation or a change in strength of the poles. Over the course of 
the simulation, the Voronoi diagram appears to retain the history of the flow held 
and deforms in the wake of the vortex pair. We also show that the relative energy 
and potential vorticity error does not exhibit secular drift. 

Before describing each of the experiments in more detail, we mention a few per- 
tinent details common to all experiments. The smoothing parameter a is choosen 
to be at least equal to the reference grid scale dx- that is the grid width of the 
initial uniform configuration of the particles. For the ID simulations, the initial 
conditions were sufficiently smooth and the approximation sufficiently high in res- 
olution that the effect of regularization was marginal and not required to ensure 
stability. However, for the 2D simulation in Experiment 3, computational cost is 
a consideration and the simulation was run at too coarse a resolution to remain 
stable unless a was increased to quadruple the reference grid size. 

7.1. Experiment 1: conservative properties of VFL. Unless otherwise stated, 
we use 128 cells, where each cell is initially of the same size. The parameters of the 
simulations are provided in Tabled] below. 





Parameter 


Value 






Number of cells N 
Time step 
Domain length 
Initial conditions 

fo 

9 

a 


128 
0.01 

2tt 

m„(t ) = 1 + 0.1exp(-800(X Q (t ) - L/2) 2 /L 2 ) 

U a {t )=0 

V a (t ) = 

2tt 

4tt 2 

dx 





Table 1. This Table lists the simulation parameters for experi- 
ment 1, which investigates the conservative properties of the VFL 
method for ID rotating shallow water. 



Results. The following Figures show the relative energy, total potential vorticity 
(PV) and total potential enstrophy PE (sum of the square of the cellwise potential 
vorticity) errors. Figure [2] shows the scaling of the energy, PV and total vorticity 
error with the size of the time step. All three quantities montonically decrease as 
the time step decreases. The error in the energy scales with the order of integrator, 
0(At 2 ). where as the potential and total vorticity do not scale with this order. This 
suggests that the VFL method introduces error in the discrete vorticity. Figure [3] 
shows (from top to bottom) the relative energy, potential vorticity and potential 
enstrophy error over 10 6 time steps, with a time step of 0.01. This number of time 
steps is approximately equal to the number of time steps required to conduct a 100 
year climate change simulation. Each error does not drift which suggests that the 
VFL method is long-time stable and conservative. 
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Figure 2 . This Figure shows the scaling laws of the mean energy 
error, potential vorticity and total vorticity with the size of the time 
step over 10 6 time steps. The energy error scales to an order At 2 
which is consistent with the order of the Stormer Verlet integrator. 
Potential vorticity (PV) and total vorticity (TV) error do not scale 
with the order of the integrator suggesting that the order of the 
error of the PV and TV are not governed by the integrator but 
by the VFL method itself. The simulation parameters are given in 
Table [TJ 

7.2. Experiment 2: Geostrophic adjustment. This experiment uses the VFL 
method to simulate the mechanism of geostrophic adjustment of shallow water in 
a plane in which the layer depth and velocity are independent of the y-direction 
(meridional direction) 

(31) U = -gV x h + f V, V = -f a U, X = U, 

Given a geostrophically balanced layer depth H°, meridional velocity V° = 
■f Vx/i and horizontal velocity U° — 0, we perturb the layer depth h = H° + h' 
by a smooth, domain centred, gaussian function hi with support L' less than the 
Rossby deformation radius Lp. The dynamics, governed by equations (|3"TT) , arc 
observed to exhibit gravity waves which propagate energy and momentum away 
from the source leaving behind geostrophically balanced flow. 

We verify in this numerical experiment that this mechanism only occurs if the 
scale of layer depth perturbation is smaller than the Rossby deformation radius by 
performing two simulations, one in which the scale of perturbation is larger and 
one is which it is smaller than the Rossby deformation radius. The parameters for 
this simulations are provided in Table [5] 

Results. Figures H] and [S] compare two different rotating shallow water regimes dis- 
tinguished by the scale of perturbation of the layer depth. In the first case, the layer 
depth is perturbed on the scale of the Rossby radius L R = 1 and in the second, 
on a scale smaller than the Rossby radius. In the latter case, the sequence of layer 
depth graphs (shown at increasing simulation times) in Figure [5] show the gravity 
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Parameter 


Value 


Number of cells N 


512 




Time step 


0.01 




Domain length 


2tt 




Initial conditions 


m a (t ) 


= 1 + Ei-i 0.01exp(-ft(X a (i ) - L/2) 2 /L 2 ) 




U a {t ) 


= 




V a (t ) 


= -^grad( J ff°) 


fo 


2tt 


9 


4ir 2 




H 


1 




L D 


2tt 




a 


dx 





Table 2. This Table lists the simulation parameters for experi- 
ment 2 which models geostrophic adjustment when (3\ = 10(L' << 
Ld) and /?2 = 1000(X' >> Ld)- H® denotes the discrete steady 
state profile, from which the meridional velocity is initialized. 



waves which propagate from the source. The bottom profile of each of these graphs 
in this Figure also shows that a geostrophically balanced region is recovered in the 
region of the source. 

Figure [6] show the corresponding energy and potential vorticity errors over 2 
rotation units (i.e. 2 days). We firstly observe that the profiles exhibit high fre- 
quency oscillations arising from gravity waves but do not exhibit drift. There are 
also jumps in the profiles with an approximate period of 0.3 days. These jumps 
occur when the gravity waves collide (recall that the domain is periodic). Figure 
[7] shows how the discrete approximation of the layer depth, over the central region 
of the domain, converges to the steady state after 0.4 days for various number of 
grid points. The outer regions of the domain exhibit gravity waves which have 
propagated away from the source. 
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Figure 3. These graphs show (from top to bottom) the relative 
energy, potential vorticity and potential enstrophy error over 10 6 
time steps. Each error does not drift which suggests that the VFL 
method is long-time stable and conservative. The simulation pa- 
rameters are given in Table [1] 
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time: 0.27 days 

0.01 | 1 1 1 r 




Figure 4. These graphs show a snapshot of the shallow water 
taken at time t = 0.27 "days" (1 "day" is a rotational unit) for 
which the initial scale of perturbation L' » Ld, where Ld is 
the Rossby deformation radius. The top two graphs show the 
layer depth perturbation and layer depth respectively. The third 
and fourth graphs show the horizontal and meridional velocities. 
The bottom graph shows the difference of the meridional velocity 
with the layer depth gradient. Note that the source region is not 
restored to a geostrophically balanced state because geostrophic 
adjustment does not occur. The simulation parameters for this 
experiment are given in Table [2] in which ft — 1000. 
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Figure 5. These graphs show snapshots of the shallow water 
taken at times t = {0.2,0.5} "days" (1 "day" is a rotational unit) 
for which the initial scale of perturbation L' << Ld, where Ld 
is the Rossby deformation radius. The top graph show the layer 
depth perturbation and layer depth respectively. The third and 
fourth graphs show the horizontal and meridional velocities. The 
bottom graph shows the difference of the meridional velocity with 
the layer depth gradient. This difference is zero when the flow is 
geostrophically balanced. Note that the source region does restore 
to a geostrophically balanced state because gravity waves prop- 
agate away energy and momentum to restore the balance. The 
simulation parameters for this experiment are given in Table [I] in 
which j3 = 50. 
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0.5 1 1.5 2 

time (days) 

Figure 6. These graphs show the relative energy (top) and po- 
tential vorticity (bottom) error profiles of shallow water over 2 
"days" (1 "day" is a rotational unit) for which the initial scale of 
perturbation is U << Ljj, where Ljj is the Rossby deformation ra- 
dius. Both profiles exhibit high frequency oscillations arising from 
gravity waves but do not exhibit drift. There are also jumps in 
the profiles with an approximate period of 0.3 days. These jumps 
occur when the gravity waves collide (recall that the domain is pe- 
riodic). The simulation parameters for this experiment are given 
in Table U in which (5 = 50. 
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Figure 7. This graph shows the layer depth at 0.4 days for various 
number of grid points. In the central region of the domain, the 
layer depth has reached a steady state (geostrophic balance). The 
graph shows how the discrete approximation of the layer depth 
converges to the steady state. The outer regions of the domain 
exhibit gravity waves which have propagated away from the source. 
This outer region is not in geostrophic balance after 0.4 days. 
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7.3. Experiment 3: Vortex pair in 2D rotating shallow water flow. This 
experiment shows the motion of a vortex pair in a purely rotational shallow water 
flow over a doubly periodic domain without bottom topography. The flow field 
is initially smooth and geostrophic balanced - the layer thickness and horizontal 
velocity are radially symmetric about the center of the domain and the meridional 
velocity is anti-symmetric about the lines x — L/2 and y — L/2. The simulation 
parameters are given in Table |3l 

A contour plot of the regularized layer thickness is shown in Figure [5] in which 
each snapshot is taken at increments of 40 days and is viewed from left to right, 
starting at the top of the page. The vortex pair is observed to undergo pure rotation 
without significant shape deformation or a change in strength of the poles. The 
Voronoi diagram corresponding to these layer thickness snapshots is shown in Figure 
[9] The diagram retains history of the flow field and largely deforms in the wake 
of the vortex pair to exhibit some intriguiging rotational structure with loss of the 
radial symmetry. Over the course of the simulation, an increasing number of cells 
deform from their initial state and the resulting effect is a more homogeneously 
deformed diagram. 

Figure [10] shows the relative energy and potential vorticity error over a period 
of 1 x 10 5 time steps or approximately 50 years. The relative energy error does 
not exhibit secular drift and is an 0(At 2 ). The relative potential vorticity error is 
within 5% and does not exhibit secular drift either. 



Parameter 


value 


N t 


1 x 10 5 (approx. 50 years) 


CO 


10 


fo 


2tt 


dt 


5 x 10~ 2 


L 


2tt 


H 


1 


N 


64 x 64 


n x 


32 


n y 


32 


a 


Adx 



Table 3. Experiment 3: the simulation parameters for the rotat- 
ing shallow water equations over a doubly-periodic domain initial- 
ized by perturbing a geostrophically balanced vortex pair. 
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Figure 8. This graph shows the layer depth at increments of 40 
days during simulation of a vortex pair in rotating shallow water 
over a doubly periodic domain. The simulation parameters for this 
experiment are given in Table [31 
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Figure 9. This graph shows the Voronoi diagram at increments 
of 40 days during simulation of a vortex pair in rotating shallow 
water over a doubly periodic domain. The simulation parameters 
for this experiment are given in Table [3] 
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Figure 10. (top) The relative energy error and (bottom) the rela- 
tive potential vorticity error over 1 x 10 time steps (approximately 
50 years). The simulation parameters for this experiment are given 
in Table GO 

8. Conclusion 

The development of new numerical methods for climate modeling is challenged by 
the simultaneous need to tractably simulate long-time dynamics and exhibit consis- 
tent potential vorticity dynamics. This paper considers the conservative properties 
of a variational formulation of the free-Lagrange method for the regularized ro- 
tating shallow water equations. The VFL method is based on the free-Lagrange 
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method, a Lagrangian method which remedies the diffusive transport computa- 
tions of Eulcrian methods, by building the transport into the evolution of the cell 
positions. Our implementation of the method has the distinguishing feature of us- 
ing both a Voronoi diagram to formulate a local mass conservation law and the 
dual Delaunay triangulation to approximate a dispersive regularization of the layer 
thickness. The semi-discrete regularized shallow water equations preserve symplec- 
tic structure and, when integrated with a symplectic time stepper, conserve energy 
over long-time simulations to the order of the symplectic integrator. 

The discrete variational principle provides a unifying point of conferred mathe- 
matical understanding from which to systematically derive new geometric methods. 
Nocthcr's theorem, for example, informs us that exact potential vorticity conserva- 
tion is only permissable when the discrete variational principle is invariant under 
continuous particle relabelling. In the absence of this property, we turn to the 
diagnostic semi-discrete potential vorticity dynamics and choose a discrete curl op- 
erator which has the property that it operates on the layer thickness gradient to 
give a zero vector field. This property guarantees the existence of a semi-discrete 
potential vorticity law. The semi-discrete potential vorticity and the semi-discrete 
divergence equations kinematically equate to the respective solenoidal and irrota- 
tional components of the velocity field. Together, they augment the description 
of how the momentum evolves in discrete space. Like the continuum divergence 
equation (expressed in the Lagrangian frame), the semi-discrete divergence equa- 
tion exhibits a div 2 U term which indicates that the flow has a very strong tendency 
to reach a purely rotational state and therefore has a marked effect upon potential 
vorticity evolution. 

Numerical results show (i) the preservation of shape and strength of an initially 
radially symmetric vortex pair in purely rotational regularized 2D shallow water 
and (ii) how the Voronoi diagram retains the history of the flow field and (iii) 
that energy is conserved to C(A 2 ) and potential vorticity error to within 5% with 
no clearly deciperable secular growth over a 50 year period. In the ID case, we 
show that the potential vorticity and energy error monotonically decreases with a 
decreasing time step, the latter as a C(A 2 ) scaling law. Simulation of a geostrophic 
adjustment mechanism provides further evidence that the VFL method exhibits 
consistent potential vorticity dynamics. 
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A.l. Calculation of the semi-discrete vorticity equations. 

(32) 

k • curl(U a ) 
= — Vt ft -dr ft 

A Z_^l a a 

^ 

- i (i XX • **) + i§w X>2- • *i - £ E i* • 

- 5 XX • + 3jS (J °»E u « ' ** + i E "f • ^ 

s v ' 

=0 

= (^Hixxx 

= R. (k • curl(Ua)) + div(U Q )k • curl(Ua). 



Remark A. 1.1. The vanishing under-braced term in the fourth line follows from 
the definition of AU^ and the property that the terms cancel around a closed loop 



(33) £u£(U^-ie-)=0- 



Note that in the continuum limit, the expression takes the form 



(34) n//"!'- - 
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A. 2. Calculation of the semi-discrete divergence equations. 

(35) 

div(U a ) = -J-£u£.dn^ 

i 

- s (£ £ ** • dn - ) + if £ <' 4 *> S u » ■ d "» - i E °* • s d "» 

- h (i ?** • -*) + E <* ■ - i ? • A < u± » 

S v ' 

-(s + ^)i?'«-*4- r - 
-(5 + i^)i? D f-**- r - 

K^ +div(Ua) )i? u «' dn ^ r - 
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